Damaging and Cracks in Thin Mud Layers. 
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We present a detailed study of a two-dimensional minimal lattice model for the description of mud 
cracking in the limit of extremely thin layers. In this model each bond of the lattice is assigned to a 
(quenched) breaking threshold. Fractures proceed through the selection of the part of the material 
with the smallest breaking threshold. A local damaging rule is also implemented, by using two 
different types of weakening of the neighboring sites, corresponding to different physical situations. 
Some analytical results are derived through a probabilistic approach known as Run Time Statistics. 
In particular, we find that the total time to break down the sample grows with the dimension L of the 
lattice as even though the percolating cluster has a non trivial fractal dimension. Furthermore, 
a formula for the mean weakening in time of the whole sample is obtained. 
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I. INTRODUCTION 

In this paper we present a careful and detailed study of 
a minimal fracture model that has been introduced at the 
aim of describing the main features of paint dessication- 
like phenomena j^] . The purpose of this work is to focus 
on the statistical properties of these phenomena on the 
basis of a recent experimental work ||^. Following the 
results of this work, we assumed that the main source of 
stress is given by the local friction between the layer of 
material and the bottom surface of the container. More- 
over, it has been noticed that the characteristic size of 
crack patterns varies linearly with the layer thickness. In 
the limit of zero thickness crack patterns lose their polyg- 
onal structure (the characteristic size of the polygons is 
zero) and become branched fractals. 

In order to model this behavior, a minimal automa- 
ton lattice model, inspired by Invasion Percolation 
and by the vectorial and scalar models described in 
has been recently introduced by the authors [Q. Here, 
we present an extended report of the study described in 
0, with a detailed description of the analytical calcu- 
lations a new numerical and theoretical results. All the 
models for quasi-static fractures, describe crack evolution 
through a non-local Laplacian field (electric field, electric 
current) acting on a solid network of bonds or sites 
In others the stress field evolves by keeping minimum the 
energy of the system. In such a case the components of 
the obtained vectorial equations are similar to the equa- 
tions describing the action of a Laplacian field jj]. In 
this model, instead, no explicit field is present. The effect 
of the stress is represented by an extremal breaking rule 
and local random breaking thresholds: at each time step, 
the bond with the smallest threshold is removed from the 
lattice. Short ranged correlations are introduced through 
a damaging of the non-broken nearest neighbors bonds 
of the just removed bond. According to the kind of frac- 
ture we deal with, one can introduce different types of 



damaging. In this paper two limiting cases are studied. 
This model is inspired by the above cited experimental 
observations |^ that, in an extremely thin layer of mud 
or paint, the only source of stress is the local friction with 
the container. Moreover, since the drying mud is a mix- 
ture of a liquid and a solid (usually amorphous) phase, 
no long range stress relaxation is present, although the 
growing crack can affect the properties of the medium in 
its neighborhood. Some important physical properties of 
the model are explained by using an approach based on 
the Run Time Statistics (RTS) scheme In particular, 
we are able to compute some relevant quantities, such 
as the evolution of the breaking probability, and of the 
probability distribution of breaking thresholds. 

The paper is so organized. In Sec. II the model is de- 
scribed. In Sect. Ill the results of numerical simulations 
are presented. In Sect. IV the model is studied analyti- 
cally and theoretical and numerical results are compared. 



II. THE MODEL 

A square lattice is considered and a quenched random 
variable Xi is assigned to each bond i. The Xj's are in- 
dependently extracted from an uniform probability den- 
sity between and I. At each time-step t, the unbroken 
bond in the lattice with the lowest value of the variable 
is broken (removed). Then damage (weakening) is ap- 
plied, in a way explained below, to the unbroken nearest 
neighbors (n.n.) of the just removed bond. After hav- 
ing introduced the damaging, the breaking and damaging 
steps are repeated until a connected, percolating, subset 
(infinite cluster) of removed bonds appears, dividing the 
system into two disconnected parts. 

Before explaining the definition of damaging, it is nec- 
essary to introduce some notations. The set of broken 
bonds up to time t is indicated with Cj, and the set of 
non-broken bonds with dCt- The number of bonds be- 
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longing to Ct is \\Ct\\ = t, while \\dCt\\ =N-t (where N 
is the total number of bonds in the lattice), in fact dCt 
is composed by the whole lattice minus the bonds in Ct- 
The definition of Ct is independent of the model. That of 
dCt, instead, can differ a lot from a model to another; for 
instance in Invasion Percolation (IP) it is simply given by 
the set of nearest neighbors of the bonds in Ct ■ 

Two different kinds of weakening of the unbroken 
neighbors are studied: Either by direct weakening or by 
re-distribution of the "stress". In the first case (rule 1), 
the unbroken n.n. are weakened, by extracting a new 
threshold x'^ between and the former value Xi. In this 
case an average weakening of one half of the former value 
at time is obtained. In the second case (rule 2), instead, 
each neighbor has a threshold weakened by a fraction 
of the threshold of the bond just removed. Both cases 
mimic the damaging produced by the enhancement of 
the stress nearby crack tips: the first case refers to a sit- 
uation where stochasticity (thermal fluctuations) is im- 
portant in the determination of the new thresholds 0, 
the second case refers to a deterministic effect around 
the crack tip. From the point of view of mud cracking, 
the two-dimensional lattice represents a very thin layer 
of mud (or paint), and quenched disorder accounts for 
local stress induced by inhomogeneous desiccation of the 
sample. Since the evolution of cracks in mud desicca- 
tion is assumed to be a slow process, the dynamics is 
assumed to be quasi-static, i.e. one microscopical break- 
ing with relative damaging for each time-step. Some au- 
thors correctly point out that otherwise time-dependent 
effects and a non-equilibrium dynamics are relevant in 
crack propagation Q . 

In this model the explicit presence of an external field 
(applied stress) and of the response of the material (strain 
of bonds) have been eliminated. The only quantity 
present is the breaking threshold, the dynamics of which 
is chosen to reproduce the evolution of cracks. This simu- 
lates the presence of a local stress field, acting not on the 
boundaries but directly on each bond. Our assumption 
is based on the experimental results in Ref. where, as 
the mud layer becomes thinner, only the inhomogeneities 
drive the nucleation of cracks. Furthermore, the hypoth- 
esis of crack developing under the same state of strain not 
only is usually applied in the presence of thermal gradi- 
ents |9|, but is also commonly reported in experiments 
of loading of softened material [l0|-[l^]. Hence, such a 
model is particularly suitable to describe, for example, 
paint drying, where the stress applied to the painted sur- 
face depends on the local action of external conditions 
(density gradient in the paint). Moreover, its simplicity 
allows us to study analytically its properties, which is a 
non common feature for fracture models. 



III. NUMERICAL RESULTS 

Numerical simulations, with cylindrical symmetry (pe- 
riodic boundary conditions in the horizontal direction) 
for various system sizes L have been performed. The 
dynamics stops as soon as a crack spanning the system 
in the vertical direction appears. Both damaging rules 
are implemented, and they are discussed in parallel. De- 
spite the simplicity of the dynamical rules, the results 
are rather interesting. We have computed the fractal 
dimension of the percolating cluster, the distribution of 
the size of clusters of broken bonds, the avalanche size- 
distribution (in order to check if long range temporal 
correlations are present), and the probability distribution 
of the breaking thresholds at the percolation time. An 
avalanche can be defined as an ensemble of causally and 
geometrically connected breakdowns (see below for a rig- 
orous definition) . Under this respect the size-distribution 
of such avalanches represents the probability of a large 
or small response of the system to an external solicita- 
tion. For example a power law distribution represents a 
critical state of the system where the response has not a 
characteristic size. 

The fractal dimension Dfoi the percolating cluster is 
computed using the box-counting method. The analysis 
is restricted to the spanning cluster to reduce the finite 
size effects present for the smaller clusters. The results 
of the box-counting analysis are reported in Table | for 
the different sizes and for the two damaging rules. The 
values oi D f for the two damaging rules coincide within 
the error bars. 

The connected clusters of broken bonds are identified 
with a standard cluster counting procedure, based on the 
Hoshen-Kopelman algorithm Also the distribution 
of finite clusters is nontrivial, showing a clear power law 
with exponent Tc = 1.54(2) (see Fig. 0(a)) for rule 1 and 
Tc = 1.57(3) for rule 2. The plots labelled with (b) in 
Fig. |l] refer to the avalanche size distribution. This quan- 
tity is interesting with respect to recent experiments p^ ] 
and models where a power law behavior of the acous- 
tic emission has been related to Self-Organized Criticality 
(SOC) [||. The presence of a SOC-like behavior would 
mean that the dynamics of fractures itself leads the sam- 
ple to a steady state where small variation of the external 
field can trigger reaction at any length-scale. In partic- 
ular the external field in this case is the applied stress, 
and the response of the sample can be considered as the 
energy released (acoustic emission) by one avalanche of 
cracks, where avalanche means a causally and geometri- 
cally connected series of breakdowns. In this oversimpli- 
fied model the external stress can be considered constant, 
since the only change after any single breakdown is due 
the damaging of the n.n.. Consequently, in this work the 
size of an avalanche is monitored as a measure for the 
acoustic emission. An avalanche can be defined as fol- 
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lows. Let us suppose that a bond i grows (i.e. it is bro- 
ken) at time t; this is the initiator of an avalanche, which 
is defined as the set of events geometrically and causally 
connected to the initial one (bond i). "Causal" connec- 
tion refers to the weakening following any bond breaking. 
In particular, when bond i grows at time t, the avalanche 
goes on at time i -I- 1 if a unbroken first neighbor bond 
j of i is removed. At time t + 2 the avalanche goes on if 
a bond k grows where fc is a unbroken first neighbor of i 
or j and so on. A linear-log plot of the probability dis- 
tribution of avalanche size, versus sample size L is shown 
m Fig. |l|b. After a power law transient, an exponential 
distribution is reached, indicating that a characteristic 
size exists for the avalanches. One can note that both 
for weakening rule 1 and weakening rule 2 simulations 
give qualitatively similar results, although for rule 2 the 
characteristic time of avalanches is smaller. This is easily 
explained, since the damaging rule 2 is less strong than 
rule 1, and, consequently, the causal connection between 
subsequent breaking events is weaker. 

This result for avalanches is similar to those obtained 
for a scalar model of Dielectric Breakdown, but differs 
from the avalanche behavior in models of fracture [§,1). 
The explanation of this behavior is motivated by two ar- 
guments. Firstly, in the present definition of an avalanche 
the threshold is changed only for the n.n.. This intro- 
duces a typical length scale, while other definitions con- 
sider as the threshold the ratio between local field and 
resistivity, thus giving the possibility of large scale cor- 
relations. Secondly, in this model broken bonds are re- 
moved from the system. This represents a substantial 
difference with many SOC models with quenched disor- 
der presented in the literature. For example, in a simple 
toy model of SOC due to Bak and Sneppen (where 
a similar refresh of thresholds is present) the dynamics 
produces clear power laws in the avalanche distribution. 
There, each site (species) deleted is replaced by a new 
one and is not definitively removed. In our model, in- 
stead, the number of candidates dCt to be broken at each 
time-step decreases in time. This is a crucial point, since 
indeed power law behavior in the presence of a scalar 
field seems to be related to a "reconstructing rule" that 
allows one to deal with a system where removed bonds 
are replaced by new ones. Therefore, only in the case 
of plastic deformation, one is in the presence of a steady 
state, as correctly pointed out by Ref. 

Therefore, the fractal dimension, the cluster size distri- 
bution and, to some extent, the avalanche size distribu- 
tion seem to be universal with respect to the two different 
local damaging rules. In the next section the study will 
focus on some quantities which, instead, are not universal 
and reproduce the evolution of the mechanical properties 
of the material during the fracturing dynamics. These 
quantities are the average probability density of break- 
ing thresholds, or histogram, (/jtix), and as a by-product 
the mean breaking threshold {x){t), which expresses the 



average resistance to breaking, or rigidity, of the system 
at time t. These quantities will be studied both numeri- 
cally, and analytically, by using a probabilistic tool called 
Rim Time Statistics (RTS) §]. 

IV. RTS DERIVATION OF THE AVERAGE 
WEAKENING OF THE MATERIAL 

As seen above, the evolution of the crack is described 
by a quasi-static extremal dynamics in a medium with 
quenched disorder. The most important question for a 
theoretical comprehension of the model is: which is the 
source of the spatio-temporal correlations developed by 
the dynamics? As pointed in |^ in relation to Invasion 
Percolation (IP), the source can be found in the memory 
effects developed by the evolution of the dynamics itself 
via an interplay between dynamical rules and quenched 
disorder. 

This can be simply realized observing that the knowl- 
edge of the growth history up to a time t, provides infor- 
mation about the probability distribution and the cor- 
relations of the random bond-thresholds. This informa- 
tion has to be added to the original information that 
the thresholds are independently extracted from the uni- 
form probability density in the interval [0, 1]. Moreover 
this information influences the probabilities of the dif- 
ferent possible continuations of the dynamics for larger 
time. This memory effects can be studied using carefully 
the notion of conditional probability. This kind of ap- 
proach to growth dynamics with quenched disorder has 
been developed in |]6|,|l7|, with particular reference to IP. 
This peculiar probabilistic algorithm is called Run Time 
Statistics (RTS). A particular modification of this tool is 
presented here taking into account the damaging mech- 
anism, which is not present in IP-like models. Finally, 
RTS is used in order to predict analytically some relevant 
quantities as the evolution of both the average probabil- 
ity density of breaking thresholds of unbroken bonds and 
of the mean resistance to breakdown x{t) of the material. 

Here we provide directly the final RTS formulas to- 
gether with a brief sketch about their meaning. A de- 
tailed derivation of the analytical results of this section is 
given in Appendix A. The RTS approach permits mainly 
to answer the following two questions, once given a cer- 
tain time-ordered geometrical path followed by the dy- 
namics up to time t: 

1. which is the effective probability density function 
of the variables Xi of the lattice conditioned to the 
knowledge of this fixed past dynamical history; 

2. which is the conditional probability of any further 
growth event at the next time-step. 

In order to introduce operative formulas, let us think to 
know the "one-bond" effective probability density func- 
tions Pi^t{x) (conditioned to the past dynamical history) 
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for each non-broken bond i. As it is clarified in Appendix 
A, this "one-bond" formulation of RTS is an approximate 
of the rigorous one. However, as shown in Q, it is a good 
approximation when the number of random numbers is 
large (as in this case). 

First of all one can write Q the breaking probability 
fii^t for each bond i at that time-step: 



dCt 

n 



dypkAy) 



(1) 



where dCt is the whole set of unbroken bonds. Note that 
at time t the number of bonds in dCt is {2L^ — t), i.e., the 
total number 2L^ of bonds in a square lattice of side L 
minus the number of broken bonds before time t. Eq. (|^) 
expresses nothing else than the effective probability that 
Xi is the minimum in the set dCt conditioned to the past 
history. The next important step is to update each pj^t{x) 
by conditioning them to this latest growth event. In this 
way one obtains the pj^t+i(x)'s conditioned to the history 
up to the time-step t+1. The effective probability density 
at time < + 1 of the latest grown bond i is usually called 
mi^t+i{x) , in order to distinguish it from the densities of 
still unbroken bonds. It is given by 



dCt 



n / dypk,t{y) 



(2) 



Eq. (||) (multiplied by dx) gives the "effective" probabil- 
ity that X < Xi < X ~\- dx, conditioned to the past fixed 
dynamical history (time-ordered path) up to time i -I- 1: 
the "memory" of the history up to time t is "recorded" 
in the set of functions Pk,tix), where k runs over all the 
bond belonging to dCt , while the last step is recorded in 
the particular functional relationship between mi^t+i{x) 
and the set {Pk,t{x)} itself. This relationship is imposed 
by the order relation among the interface variable Xk, i.e. 
by the fact that Xi is the minimum in dCt- Note that, 
once a bond is broken, it does not participate anymore 
to the dynamics. For this reason, the "effective" prob- 
ability density function of its threshold does not change 
anymore in time and is given definitely by Eq. §). 

For the remaining bonds one has to distinguish among 
the unbroken bonds far away from the bond i and the 
unbroken nearest neighbors bonds, which will be weak- 
ened by the growth of bond i. The updating rules, for 
the two different mechanisms of damaging, differ only for 
this last set of bonds. For the non- weakened bonds, one 
has in both cases the following updating equation: 



1 f'-" 
Pj,t+i{x)^ PjA^) dyp,Ay) 



dCt ^1 

W / dzpkAz) 

k(^i,3)-'y 



(3) 



The updating equations for the weakened bonds are 
instead the following: 
(1) For the damaging mechanism 1: 



Pj,t+i(x) = / dy -9{y - x)pj,t{y) x 

i^i,t Jo y 



X / dzpi,t{z) 
Jo 



dCt „i 

W j dupk^iu) 



(4) 



(2) For the damaging mechanism 2 (see Appendix A): 



Pj,t+i{x) = 

^ „1 dCt „l 

= dy TT / dzpk,t{z) 



Pi.tiy)Pj,t[ X H I X 



xe 



riif 



rii.t - 1 



x-y]0 (ni,f(l ~x)-y) 



(5) 



Note that the main difference between Eq. (||) and Eq. (|^) 
is due to the fact that the number of n.n. rii^t of the bond 
i at time t appears explicitly only in the latter, i.e. only 
in the second model (rule 2) the damaging is an explicit 
function of the geometry, while in the former (rule 1) the 
damaging is a "one-bond" process. 

Eqs. (|l|-||) coincide with the ones introduced for the 
RTS approach to IP (apart from the different definition 
of the growth interface dCt). Eqs. (^, instead, are 
new and account for the n.n. weakening. Eqs. (Qj^) allow 
one to study the extremal deterministic dynamics as a 
kind of stochastic process with memory. In particular, 
can be used to evaluate systematically the statistical 
weight of a fixed time-ordered growth path, while the 
Pj^t(a;)'s store information about the growth history. 

A very important quantity to characterize the prop- 
erties of the dynamics is the empirical distribution (or 
histogram) of unbroken thresholds. This quantity is de- 
fined as: 



(6) 



where, ht{x)dx is the number of non- broken bonds be- 
tween X and a; -|- c?a; at time t, conditioned to the past 
dynamical history. 

Considering the effect of the growth of bond i at time 
t on this quantity, one gets 



ht+i {x) = ht {x) - mi^t+i {x) Pj,t (^) + y^Pj,t+i(^') 

3{i) Ki) 



(7) 



where indicates the sum over the rii^t unbroken n.n. 
of i. Moreover, rrii^t+iix) aiid Pj,t+i{x) are given respec- 
tively by Eq. (|) and Eqs. (|), (|) ( (|) for rule 2) . 
Being the histogram an almost self-averaging quantity of 
the model in the large time limit, one can evaluate its 
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shape in the "typical" reahzation of the dynamics tak- 
ing the average over all the possible histories up to time 
t+1. The notation (...) is introduced to indicate for this 
average. The l.h.s. of Eq. (|) can be computed as 

{ht+,ix)) ^\\dCt+i\\^t+iix)^[N-it + l)]<f>t+iix), (8) 

where iV = 2L^ is the total number of bonds in the lattice 
and (x) represents the average thresholds density func- 
tion over the unbroken bonds at time t (normalized to 1), 
i.e. (f>t{x) — Pk,t{x) where fc is a generic interface bond. 
For the r.h.s. of Eq. (^ the main difficulty arises in the 

evaluation of {nii^t+i) and (j2j{i)Pj-t+ii^)^- Following 
[pi , one can write 



(m,,t+i) ~ {N -t)(l)tix) 



1 - / dy<j)tiy) 





N-t-l 



(9) 



In obtaining Eq. (|^), we used the definition of (j)tix) and 
the following approximation: 



Y[ Pk,t{xk)) = Yi {Pk,t{xk)) = Yl (l^tixk)- 
\kedCt I kedCt kedCt 



Using again the definition of 0t (x) , one gets 



(10) 



(11) 



where nt — {rii^t)- Using Eq. (^), corresponding to 
the weakening rule 1, and the approximations given by 
Eq. (p^), one has 



V N-t-l .L ^ y 














1 - 


f dz(f)t{z) 








Jo 





(12) 



The equation for the (t>t+i (x) for rule 1 will finally read: 

, , N -t-nt , , 
(IJt+i(x) = — -Mx) + 



N-t-l 



N-t 
'N-t-l 



4>t{x) 



1-1 dy(t)t{y) 



N-t-l 



N-t [\ My) ^ 

-nt , — — I ay x 



[N-t-iy 



y 









{- 


1- r dzMz) 


N-t-l"^ 




Jo 





(13) 



Note that even at percolation time iV — t is a large 
number. For this reason terms in Eq. (13) containing 



the term [l — dycfit (j/)] ^ ' ^ are negligible for x such 



that dy(j)t{y) is finite (i.e. larger than 1/{N — t)). It 
is easy to show that the continuum limit of Eq. (|l^) , for 
such values of x, is invariant under the rescaling L aL 
(i.e. TV — > a^N) and t — > a^t. This result is based on the 
assumption that: 



nt{L) = na2t{aL) 



(14) 



The numerical simulations suggest the following scaling 
form for nt{L) (see Fig. ph: 



nt{L) 



1 



1 + t/AL^ 



(15) 



where (3 = 0.23(2), A = 0.030(2) and n^aa; = 6 is the 
lattice coordination number. This form for nt satisfies 
Eq. (0). 

The study of the weakening rule 2 is quite similar. 
Eqs. (|l]-^) keep the same, while the conditioned proba- 
bility density for the weakened bonds is given by Eq. (||). 

By following the same steps as above, the following 
equation for the (l)t+i{x) is obtained: 



(t)t+i{x) = 
N-t 



N -t-nt 
N-t-l 



(t)t{x) 



N-t-l 
N-t 



(t)t{x) 



dyMy) 



N-t-l 



N-t-l 



-nt 



xcpt \ x-\ 

nt 



dy 



1 - 



dz(j)t{z) 



N-t-2 



My) X 



nt 



nt 



X - y ]9{nt{l - x) - y) (16) 



All the assumptions we made for the case 1, including 
the scaling ansatz given in Eq. (|l^, are valid for case 2. 
In particular, from numerical simulations, one can find 
the following behavior for nt{L) (sec Fig. ||) : 



nt{L) 



: exp 



(17) 



The analytical study of both kinds of weakening allows 
to make three important predictions: (1) Firstly, we 
find both theoretically, from the numerical solution of 
Eqs.(^, [l^ ), and from the numerical simulations of the 
model, a discontinuity in the histogram (see Fig. ^ and 
^ , indicating that the system evolves in such a way as to 
remove all bonds with threshold smaller than some crit- 
ical value Xc.{2) Second, from the symmetry properties 
of Eqs. ([l|, hq), we deduce that the number tsp{L) of 
broken bonds at the percolation time is proportional to 
i^, even though the percolating cluster is fractal. This 
result, confirmed by numerical simulations (see Fig. |qa , 
Fig. 1^), and compatible with the scaling function ( |l5| ) 
for nt{L), is deduced supposing that at the percolation 
time the shape of the histogram is independent of i, an 
assumption which fits well with the numerical histogram 
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(see Fig. ^ and Fig. ||a). (3) Finally, we present an ap- 
proximated result for the dynamical behavior of the av- 
erage value (over the unbroken bonds) of the thresholds 
(x) (t). This quantity can be seen as a characterization 
of the average resistance of the material in time. 

In order to find the evolution equation of (x) (t) for the 
damaging rule 1, it is enough to multiply both sides of, 
respectively, Eq. ([l^) and Eq. (|l^) for x and integrate 
them in the whole interval [0, 1]. Then one finds: 



l + nt/[2{N -t-1)] '■^ 



c)(t) 



N -t-1 



dx 



1 - 



-I N-t 



dy4>t{y) 



(18) 



For the damaging rule 2, the way to find the equation for 
{x) is even simpler. In fact, it is enough to consider that 
at each time step, the global effect on {x) is equivalent 
to remove two bonds with the resistance equal to the 
minimal one at that time. Therefore, one can write: 



{x)(t + l)=[l + 



1 



N -t-l 



{x) it) 



N -t-1 



dx 



n N-t 



dyMy) 



(19) 



For rule 1, it is simple to see, from Eq. (|18[), that 
{x){t + 1) < {x){t) until rif > 2 (which is verified for 
all the times). This means that on average the medium 
weakens during the evolution even if the weakest bond is 
removed at any time step. This is due to the fact that, in 
this case the weakening of the neighbors of the weakest 
interface bond has a stronger effect on the material than 
the removal of the weakest bond itself. For the rule 2, in- 
stead, one finds that {x){t + l) > {x){t), if (x) (t) is larger 
than the double of the average minimal threshold, and, 
due to the extremal nature of the dynamics, this is veri- 
fied always in the large N limit, i.e. in the limit of a large 
number of bonds in the interface at any time-step. This 
means that in this second case, damaging is not strong 
enough to allow a global weakening of the system, which 
becomes more and more rigid. This is reasonable since 
in rule 2 the stress on the weakest bond is redistributed 
to the nearest neighbors, and the total initial stress is 
conserved, while in rule 1 there is not total stress con- 
servation. In other words, in the model with rule 1 the 
damaging is a multiplicative effect, i.e. the damaging 
is proportional to the old threshold (which can be big), 
in the model with rule 2 the damaging is quite reduced 
by the fact that at each time-step it is proportional to 
the minimal threshold in the whole system. In Figs. 
^ the time evolution of {x){t) obtained from computer 
simulations is compared with the theoretical prediction. 
Our analytical results are in good agreement with nu- 
merical simulations. For rule 2, numerical simulations of 
the histogram evidentiate a low x tail below the critical 



threshold, which tends to disappear as the system size 
grows, and a non zero slope of the part just above the 
critical threshold. The first one is a clear finite size ef- 
fect, which is less important in the simulations of rule 1, 
because for rule 1 the critical threshold is very small. Of 
course, such a finite size effect is absent in the theoret- 
ical results, as in all mean field (MF) approaches. The 
second effect could be due to spatial correlation induced 
by the damaging rule 2, which in the analytical approach 
are neglected. This second effect does not disappear as 
the system size grows. Consequently the agreement be- 
tween the numerical simulations of (x) (t) and Eq. (^9|) is 
less good than for rule 1. The numerical {x)(t), mainly 
because of the nonzero negative slope of (/)(a;)above Xc, is 
a bit smaller than the theoretical prediction. 

With respect to real fracturing processes the behavior 
of the average resistance {x) {t) obtained with rule 2 is 
more realistic, since in real materials one usually observes 
that the material during micro-cracks formation becomes 
more rigid, although more fragile, since the number of 
bond one has to break to have global breakdown becomes 
smaller and smaller. Moreover it is worth to note that, 
apart the shape of (j>{x) and the behavior of {x){t), all the 
other statistical properties of the system do not depend 
on the used weakening rule. 

Finally, it is worth to point out that, to our knowledge, 
apart the qualitative results of j|] , no quantitative experi- 
mental results are available. For example, a measurement 
of the fractal dimension of cracks or their size distribution 
would be extremely useful to further test the predictions 
of this model. At the moment, this model seems able 
to capture, with its extremely simplified dynamics, some 
basic properties of fracturing processes. 

In conclusion, we have presented a new model for frac- 
tures, which is useful in describing in a semi-quantitative 
way some basic mechanism in drying paint-like and mud- 
like cracking processes, for extremely thin samples. Due 
to its extreme simplicity, the model is particularly suit- 
able for large scale simulations and takes into account the 
damaging effects involved in fracture propagation. Even 
in this simple model we are able to analyze which con- 
ditions trigger SOC behavior in such systems. Further- 
more, the change in the threshold distribution, induced 
by the damaging mechanism, allows us to write down 
explicitly the form of the breakdown probability for the 
bonds of the sample. Possible further research could con- 
sist, for damaging rule 2, in a more refined calculational 
scheme, in which two variable probability densities are 
also considered. This would be the first order correc- 
tion to our MF approach considering only one-variable 
distributions, and could allow us to take into account 
correlations induced by the damaging rule. Such a gen- 
eralization of the RTS theory, formally discussed in , 
is however technically very difficult. Another research 
direction we are following is the application of real space 
techniques, combined with the RTS approach, to calcu- 
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APPENDIX A: RTS FOR THE DAMAGED 
SYSTEM 

In this appendix a simple explication of the RTS prob- 
abilistic equations is provided. 

First of all, one has to note that in this kind of models 
(as well as in Invasion Percolation) the initial condition 
of the system is characterized by independent variables 
(the breaking thresholds of the bonds) identically and 
uniformly distributed. 

However, once the minimal value in the set is found 
and the relative bond broken, the knowledge of this 
event makes the variables of the remaining non-broken 
bonds no longer simply uniformly distributed in the in- 
terval [0,1], and correlated (no more independent one 
each other). In fact, after the breaking of the bond with 
the minimal threshold, one has to condition the proba- 
bility of any further event to the last known event. This 
information influences the probability distribution of the 
remaining bonds of the system and creates correlations 
among them ]l9| ] . 

The systematic study of this "memory" effect is what 
is called Run Time Statistics |p|,p7t. 

In order to clarify the "step by step" mechanism of 
storage of conditional information, let us think to have 
fixed a time-order path At, i.e. an history of the dynamics 
up to time t. At is given by the time ordered sequence 
{io, ?2, «t-i} of the broken bonds up to time t. Let us 
suppose to know the joint threshold probability density 
function Pt{{x}gct\At) of the whole set of non-broken 
bonds conditioned to the knowledge of the past history 
At- Pt{{x}dCt\^t) represents the "effective" distribution 
of the disorder at the t*^ time-step of a fixed history At . 
Note that at time t ^ 0. one has: 



Po{{x}dCo) ^Yipoixk) ^1, 



(Al) 



k£S 



where S is the whole lattice, as no information from the 
dynamics is still present. 

Since any kind of "order" relation, superimposed to a 
set of independent stochastic variables, introduces cor- 
relations, in general Pt{{x}QCt\-^t) does not factorize in 
the product of single-bond "effective" density functions 
for t > That is, it is not possible to write: 

Pt{{x}acMt) ^ II PkA^k) . (A2) 
kedCt 

However, as shown in [ p^ , in the limit of large number of 
variables the "geometrical" correlations m Pt{{x}gct\^t) 
become negligible, and one can make, at any time step. 



the approximation given by Eq. ( A2 ) . Therefore, we con- 
sider the approximated case where the "effective" proba- 
bility density function of the disorder of the system, with 
all the information about the past history stored, is given 
by the set of "effective" one-bond functions Pk,t{x) for 
each non-broken bond k. The rigorous exposition of RTS, 
by using the non-factorizable function Pt{{x}QCt\-^t) at 
each t is given in [ pTf . 

Knowing the set of functions pk^t{x), one can write the 
"effective" probability /ii_t that a given bond i of the set 
is broken at that time. It is simply the probability, condi- 
tioned to the whole past history, that Xi is the minimum 
in the set of non-broken bond variables. Consequently, it 
is given by Eq. (|l]), i.e. 



tJ^i,t 



dxpi^t{x) 



dCt ^1 

n / d,yPk.t{y) 



(A3) 



The set of /i^^t for each non-broken bond and for each 
time-step, defines a branching process of the dynamics; 
i.e. each history At at time t continues with a certain 
probability /i^.f in a different history At+i at time t + 1 
for each breaking bond i at time t. In order to continue 
the probabilistic description of the branching at further 
time-steps, one should obtain the new set of functions 
Pk^t+i{x) for these different cases of breaking at time t, 
using only the "old" set of Pk,t{x) and the set of prob- 
abilities ^i^t defining the branching. This is possible by 
using the notions of conditional probability. Here the 
simple rule relating the conditional to joint probability 
of a first event A to a second event B is reminded p9[ |: 



Prob(A|B) = 



Prob(A n B) 
Prob(B) 



(A4) 



where, as usual, A\B means the event A conditioned to 
the event B, while Af^ B the event A joint to the event 
B. 

Note that "memory" up to time t, for a fixed history 
At in the branching of all the possible histories, is already 
stored in the functions Pk,t{x). Consequently, in order to 
obtain the set of probability functions Pk.t+i{x) for the 
history At+i obtained from At adding the breaking of 
bond i at time t, one has to store only information about 
the last step. 

At this point one has to distinguish the three cases: 
(1) the just broken bond i, (2) a non-broken bond j far 
from i, and (3) a non-broken neighbor / of i. 
(1) In this case let us call the conditioned probability 
density of bond i after its braking with nii^t+iix) instead 
of Pi^t+i{x), remarking with this that after its break- 
down, bond i is removed definitely from the interface. 
Note that, since after t the bond i does not participate 
to the dynamics, its "effective" probability density will 
not change anymore. mi^t{x)dx is the probability that 
X < Xi < X + dx, conditioned to the the past history up 
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to its breaking. However, since the memory up to the 
time-step just before its breaking is stored in the known 
functions pfe^t (a;), rai_t{x)dx is the probabihty, calculated 
using the set of functions{pfc^i(a;)}, that x < Xi < x + dx 
(event A of Eq. (A4)) conditioned to the fact that the 
bond Xi is the minimum in the set of interface bonds 
at that time (event B of Eq. (A4)). Therefore, from 
Eq. (0), one has Eq. (|): 



m^^t+i{x) = —pi^t{.x) 



dCt 1 

n / dypk.tiy) 



(A5) 



In a quite similar way we can update the effective prob- 
ability densities for the cases (2) and (3). In the case 
(2), using the set of functions {pfe^t(x)}, pj t+i{x)dx is 
the probability that x < Xj < x + dx (event A) condi- 
tioned to the fact that Xi was the minimal in the interface 



at time t. Again from Eq. (A4) one has Eq 



Pj,t+i{x) ^ PjAx) dyp^^tiy) 

fJ'i,t Jo 



W / dzpk,t{z) 



(A6) 



In the case (3) one has to distinguish the two different 
damaging rules, and the conditioning events are more 
complex. For rule 1, using the set of function {pfe j(a;)}, 
Pi,t+i{x)dx is the probability that x < xi < x + dx (event 
A) conditioned to the fact that Xi was the minimum and 
that the value of xi at this time-step differs from the 
value at the previous time-step for a random fraction of 
itself (event B). One, then, gets Eq. (^): 

1 I 

Pj,t+i{x) = / dy -6{y - x)pjAy) x 



X / dzpi^z) 
Jo 



dCt ^1 

Y\ / dupk,t{u) 



(A7) 



Finally, for rule 2, always using the set of functions 
{Pk,t{x)}, pi^t+i{x)dx is the probability that x < xi < 
x + dx (event A) conditioned to the fact that Xi was the 
minimum and that the value of xi at this time-step dif- 
fers from the value at the previous time-step for a fraction 
l/rii^t of Xi (event B). From this one has eq. (E): 



1 



dy 



^^i,t Jo 



z) 



. rii.t - 1 



dCt ^1 

n / 

x~y]0 {n.i^{l -x) -y) 



dCt ^1 

n / dzpk,t{ 



Pi,tiy)pj,t[ X + I X 



(A8) 
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L = 64 


L = 128 


L = 256 


Df (dam. rule 1) 


1.75(2) 


1.74(2) 


1.74(2) 


Df (dam. rule 2) 


1.73(2) 


1.75(2) 


1.76(2) 



TABLE L Fractal dimension of the spanning cluster for 
different sizes and for the two damaging rules. 
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• simulations, L=32 
n(t)=n,„„^exp(-t/t„(L)), t„(L)~ l' 
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FIG. 1. (a): Probability distribution {logio-logio plot) 
Dc{s) of the cluster size, for L = 64, 128, 256. (b): Avalanche 
size distribution (linear-Zogio plot) Davit) for L — 128,256. 
(c) and (d): The same quantities for the weakening rule 2. 



• simulations, L=32 
- n(t)=n^JVa+t/AL')t, p=0.23{2) 
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• simulations, L=64 



n(t)=n_[l/(l+tyAL )] , (5=0.22(2) 
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FIG. 2. fit of nt(L) with the scaling form |l^ (weakening 
rule 1) for L = 32 (a) and L = 64 (b). 



• simulations, L=64 
- n(t)=n„^exp(-t/t„(L)), t„(L)--L' 



(b) 



FIG. 3. fit of nt{L) with the scaling form |l^ (weakening 
rule 2) for L = 32 (a) and L = 64 (b). 
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FIG. 4. Solution of Eq. ^ for the histogram (t>t{x) at the 
spanning time (b), compared with simulations (a), for the 
weakening rule 1. 
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FIG. 5. Solution of Eq. ^ for the histogram (t>t(x) at the 
spanning time (b), compared with simulations (a), for the 
weakening rule 2. 
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• numerical simulations L=128 
- theory, L=128 
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FIG. 7. (a) Spanning time versus system size L for weak- 
ening 2. Also in this case, there is very good agreement with 
the expected scaling law tsp{L) oc . (b) Solution of Eq. |l^ 
compared with numerical simulations. 
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FIG. 6. (a) Spanning time versus system size L for weak- 
ening rule 1. One can see a nice agreement with the expected 
scaling law tsp{L) oc L^. (b) Solution of Eq. |l^ compared 
with numerical simulations. 
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